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ABSTRACT 

Short-period exoplanets can have dayside surface temperatures surpassing 2000 K, hot enough 
to vaporize rock and drive a thermal wind. Small enough planets evaporate completely. We 
construct a radiative-hydrodynamic model of atmospheric escape from strongly irradiated, 
low-mass rocky planets, accounting for dust-gas energy exchange in the wind. Rocky planets 
with masses < O.IM® (less than twice the mass of Mercury) and surface temperatures > 
2000 K are found to disintegrate entirely in < 10 Gyr When our model is applied to Kepler 
planet candidate KIC 12557548b — which is believed to be a rocky body evaporating at a 
rate of M > 0.1 Mgj/Gyr — our model yields a present-day planet mass of < 0.02 or 
less than about twice the mass of the Moon. Mass loss rates depend so strongly on planet 
mass that bodies can reside on close-in orbits for Gyrs with initial masses comparable to or 
less than that of Mercury, before entering a final short-lived phase of catastrophic mass loss 
(which KIC 12557548b has entered). Because this catastrophic stage lasts only up to a few 
percent of the planet's life, we estimate that for every object like KIC 12557548b, there should 
be 10-100 close-in quiescent progenitors with sub-day periods whose hard-surface transits 
may be detectable by Kepler if the progenitors are as large as their maximal. Mercury-like 
sizes (alternatively, the progenitors could be smaller and more numerous). According to our 
calculations, KIC 12557548b may have lost ~70% of its formation mass; today we may be 
observing its naked iron core. 
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1 INTRODUCTION 

Atmospheric escape shapes the faces of planets. Mechanisms for 
mass loss vary across the Solar System. Planets lacking mag- 
netic fields can be stripped of their atmospheres by the Solar 
wind: the atmospheres of both Mercury and Mars are continu- 
ously erode d and refreshed by Solar wind ions jPotter & Morgar 



Killen et al 



1990, 1997: Lammer & B auer ^ 19971 : Isida et al1l200(t 
2004. : Jakoskv & Phillips 20011) . Venus demonstrates the extreme 



sensitivity of atmospheres to stellar insolation. Although its dis- 
tance to the Sun is only 30% less than that of the Earth, Venus 
has lost its water to hydrodynam i c escape powered by Sol ar ra- 
diation jKasting & Pollackl 1 19831 : IZahnle & Kastind Il98d ). Im- 
pacts with comets and asteroids can also purge planets of their 
atmosphe res, as has thought to have happened to some e xtent 
on Mars l lMelosh & Vickervlll989l : llakoskv & Phillips! 1200 ih and 
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some giant planet satellites ( IZahnleetal]|l992h . An introductory 
overview of atmospheri c escape in the Solar System is given by 
leading & Zahnld j2009l) . 

Smaller bodies lose their atmospheres more readily because 
of their lower surface gravities. Bodies closer to their host star also 
lose more mass because they are heated more strongly. Extras olar 
planets on close-in orbits are expected to have highly processed at- 
mospheres. In extreme cases, stellar irradiation might even evapo- 
rate planets in their entirety — comets certainly fall in this category. 

Hot Jupiters are the best studied case for atmospheric 
erosion in extrasolar planets (for theoretical models, see 



Yelle 2004 20061: Ixian et al, 20051: 


Garcia Muiiozl 20071: 


Holmstrom et alj 20081: Murrav-Clav et all 


20091: lEkenback et al.l 


2010l: TrembUn & Chiand 2013). For these 


gas giants, mass loss 



is driven by stellar ultraviolet (UV) radiation which photo-ionizes 
hydrogen and imparts enough energy per proton to overcome 
the planet's large escape velocity (~60 km/s). Mass loss for 
hot Jupiters typically occurs at a rate of M ~ lO'^-lO" g/s 
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so that ~1% of the planet mass is lost over its lifetime jYelld 
I2OO6I : iGarcia Muno3 I2OO7I : iMurrav-Clav et alj |2009|) . Although 
this erosion hardly affects a hot Jupiter's internal structure, 
the outflow is observable. Winds escaping from hot Jupiters 
have been spectroscopically detected by the Hubble Space 
Telescope (HST). During hot Jupiter transits, several absorp- 
tion lines (Hi, O i , C n, Mg n, and S i m) deepen by ~2-10% 
dVidal-Madjar et al.' '200 3L |2004 I2OO8I: ISen-Jaiffell I200I I2OO8I: 
Fo ssati et al.,2010; Leca ve lier Pes Etangs et alfcOlOlJLinskv et alJ 
20ld : iLecavelier des Etangs et al l2012h . 



Because of their lower escape velocities, lower mass super- 
Earths should have their atmosph eres more significantly sculpted 
by evaporation. iLopez et al.l ( l2012h found that hydrogen-dominated 
atmospheres of close-in super-Earths could be stripped entirely by 
UV photoevaporation (e.g., Kepler-llb). 

At even lower masses and stellocentric distances, planets 
might evaporate altogether. When dayside temperatures surpass 
~2000 K, rocks, particularly silicates and iron (and to a lesser ex- 
tent, Ti, Al, and Ca), begin to vaporize. If the planet mass is low 
enough, thermal velocities of the high-metallicity vapor could ex- 
ceed escape velocities. Under these circumstances, mass loss is not 
limited by stellar UV photons, but is powered by the incident bolo- 
metric flux. 



1.1 KIC 12557548b: A Catastrophically Evaporating Planet 

In fact, the Kepler spacecraft may have enabled the discovery of 
an evapor ating, close-in ro c ky pla net near the end of its life. As re- 
ported bv lRappaport et all Eom . the K-star KIC 12557548 dims 
every 15.7 hours by a fractional amount / that varies erratically 
from a maximum of / a 1.3% to a minimum of / < 0.2%. These 
occultations are thought to arise from dust embedded in the time- 
variable outflow emitted by an evaporating rocky planet — here- 
after KIC 1255b. Because Kepler observes in the broadband op- 
tical, the obscuring material cannot be gas which absorbs only in 
narrow lines; it must take the form of dust which absorbs and scat- 
ters efficiently in the continuum. As we will estimate shortly, the 
amount of dust lost by the planet is prodigious. 

From the ~2000 K surface of the planet is launched a ther- 
mal ("Parker-type") wind out of which dust condenses as the 
gas expands and cools. Stellar gravity, together with Coriolis and 
radiation-pressure forces, then shape the dusty outflow into a 
comet-like tail. The peculiar shape of the transit light curve of KIC 
1255b supports the interpretation that the occultations are caused 
by a dusty tail. First, the light curve evinces a flux excess just be- 
fore ingress, which could be caused by the forward scattering of 
starlight toward the observer. Second, the flux changes more grad- 
ually during egress than during ingress, as expected for a trail- 
ing tail. Both feat ures of the light curve and their explanations 
we re elucidated by Rappaport et al.l ( l2012h . and further quantified 
by iBrogi et al] ( 1201 2I) who presented a phenomenological light- 



scattering model of a dusty tailQ 



' As far as we ai'e aware, no other model is as successful at explaining the 
light curve of KIC 1255b. The fact that egress is more gradual than ingress 
argues agai nst mass transfer by Roche lobe overflow from a companion to 
the star (c.f. iLai et al.ll201ol) ; the accretion stream would lead the compan- 
ion in its orbit, not trail it as is required by o bsei'vation. Bow shock mo d- 
els of planets pass ing through stellar coronae l lVidotto et al.l201(ll201 libt 
lLlamaetal.l201lh . motivated by the ingress-egress asymmetry observed for 



Spectra of KIC 12557548 taken using the Keck Telescope re- 
veal no radial velocity variations down to a limiting accuracy of 
~100 m/s (A. Howard and G. Marcy 2012, personal communica- 
tion). A deep Keck image also does not reveal any background 
blended stars. These null results are consistent with the interpreta- 
tion that KIC 1255b is a small, catastrophically evaporating planet. 

The characteristic mass loss rate for KIC 1255b — a quantity 
that we will reference throughout this paper — can be estimated 
from the observations as follows. First compute the total mass in 
grains required to absorb/scatter a fraction / of the starlight, as- 
suming that the dust is optically thin (or marginally so). Out of the 
total area nR^ presented by the stellar face, take the dust cloud to 
cover an area A and to have a characteristic line-of-sight optical 
depth "LiKi, where is the grain mass per unit area (on the sky) 
and is the grain opacity. Then 



f- 



(1) 



Assume further that grains can be modeled as a monodispersion of 
spheres with radii s and internal bulk density pi„t so that 

3 



Kd = . 

4pinti 

The total mass in dust covering the star is then 



(2) 



(3) 



To obtain a mass loss rate, divide this mass by the orbital period 
of 15.7 hours on the grounds that transit depths change by up to an 
order of magnitude from orbit to orbit. For / = 1% — when KIC 
1255b is in a state of relatively high mass loss — the resultant mass 
loss rate in grains is 



M 



1255b,dust 



0.4 



R, 



f 



0.65Re) \ 0.01 /\ 0.5 yum/\lg/cm3 



Pir 



M© Gyr 



(4) 

This is the mass loss rate in dust only. The gas-to-dust ratio by 
mass, Tj, must be at least on the order of unity if outflowing gas is 
to cany, by E pstein drag, embedded grains out of the gravity well 
of the planet l lRappaport et al.ll2012h . Henceforth we will adopt as 
our fiducial total mass loss rate (for the relatively high value of 
/ = 1%): 



Mnssb = (1 + ?/)Mi255b.dust ~ 1 Gyr 



(5) 



The total mass loss rate tx pin,.?(l -1- rj). Given the uncertainties in 
these three factors, the total mass loss rate could take any value 
^i255b > 0.1 Mgj / Gyr (but see also our Figures l8l and [TOl which 
plac e upper bounds on mass loss rates). 

IRappaport et al.l 1 I2OI2I) claimed that the present-day mass of 
KIC 1255b is ~0.1 Mq. They estimated, to order-of-magnitude, the 
maximum planet mass that could reproduce the observed mass loss 
rate given by equation (|5}. Here we revisit the issue of planet mass 



the exoplanet WASP-12b iFossati et al]|201(]t) . do not apply to KIC 1255b 
because the WASP-12b obsei'vations pertain to nan'ow gas absoiption lines, 
whereas the Kepler observations are taken in a broadband filter At the low 
gas densities characterizing exoplanet winds, only dust can generate the 
large continuum opacit ies required by the deep transits observed by Kepler. 
IRappaport et alj j2012h reviews these and additional arguments supporting 
the dust tail mo del. A similar interpretation was proposed for the light c urve 
of beta Pictoris iLecavelier Des Etangs et al.ll997HLamers et aljl997l) . 



© 0000 RAS, MNRAS 000, 000-000 



Catastrophic Evaporation of Rocky Planets 3 



— and the evaporative lifetime it implies — with a first-principles 
solution of the hydrodynamic equations for a planetary wind, in- 
cluding an explicit treatment of dust-gas thermodynamics. We will 
obtain not only improved estimates of the present-day mass, but 
also full dynamical histories of KIC 1255b. 



1.2 Plan of This Paper 

We develop a radiative-hydrodynamic wind model to study mass 
loss of close-in rocky planets with dayside temperatures high 
enough to vaporize rock. Although our central application will be to 
understand KIC 1255b — both today and in the past — the model 
is general and can be easily be modified for parameters of other 
rocky planets. 

The paper is structured as follows. In ^we present the model 
that computes M as a function of planet mass M. We give results in 
^ including estimates for the maximum present-day mass and the 
maximum formation mass of KIC 1255b. A wide-ranging discus- 
sion of the implications of our results — including, e.g., an explana- 
tion of why mass loss in this context d oes not obey the o ften-used 
energy-limited mass loss formula (cf. iLopez et"al]l2012h . as well 
as the prospects of detecting close-in Mercuries and stripped iron 
cores — is given in ^ Our findings are summarized in ^ 




T [K] 

Figure 1. Equilibrium vapor pressures of pyroxene and olivine vs. temper- 
ature, obtained from equation 0. 



2 THE MODEL 

We construct a ID model for the thermally-driven atmospheric 
mass loss from a close-in rocky planet, assuming that it occurs in 
the form of a transonic Parker wind. The atmosphere of the planet 
consists of a high-metallicity gas created by the sublimation of sil- 
icates (and possibly iron, as we discuss in 94.7) from the planetary 
surface. We focus on the flow that streams radially toward the star 
from the substellar point on the planet. The substellar ray is of spe- 
cial interest because the mass flux it carries will be maximal insofar 
as the substellar point has the highest temperature — and thus the 
highest rate of rock sublimation — and insofar as stellar gravity 
will act most strongly to accelerate gas away from the planet. 

As the high-metallicity gas flows away from the surface of the 
planet, it will expand and cool, triggering the condensation of dust 
grains. Dust grains can affect the outflow in several ways: (a) latent 
heat released during the condensation of grains will heat the gas; (b) 
continuously heated by starlight, grains will transfer their energy 
to the gas via gas-grain collisions; and (c) absorption/scattering of 
starlight by grains will reduce the stellar flux reaching the planet, 
reducing the surface temperature and sublimation rate. In this work, 
we do not account for the microphysics of grain condensation, as 
cloud formation is difficult to calculate fr om first principles ( see, 
e.g., the series of five pap ers beginning with Helling et al .1200 ll ; see 
also lHelling et alj2008d and lHelling et alj2008bl) . Instead, we treat 
the dust-to-gas mass fraction .yjus, as a free function and explore 
how our results depend on this function. In our 1-fluid model, dust 
grains are perfectly entrained in the gas flow, an assumption we test 
a posteriori ( 94.21 1. 

We begin in g2.1l bv evaluating thermodynamic variables at 
the surface of the planet — the inner boundary of our model. In 92.21 
we introduce the ID equations of mass, momentum, and energy 
conservation. In 92.3 1 we explain the relaxation method employed 
to solve the conservation laws. 



2.1 Base Conditions of the Atmosphere 

Dayside planet temperatures peak at approximately 



T 

SI] 



T. 



-r.,„rr.cc/4 ^ 2150 g-^-'f-'--/* K, 



(6) 



where for our numerical evaluation we have used parameters ap- 
propriate to KIC 12557548 (orbital semimajor axis a x 0.013 
AU; effective stell ar temperature T, « 4450 K; and stellar ra- 
dius R, x 0.65 j^f:^; lRappaport et al.ll2012l) . Here Tsu^facc is the opti- 
cal depth between the planet's surface and the star due to absorp- 
tion/scattering by dust grains. 

We expect rocky extrasolar planets to consist of the same 
silicates predominantly found in the mantle of the Earth: pyrox- 
ene ([Fe,Mg]Si03) and olivine ([Fe,Mg]2Si04). At the estimated 
Tsmfmc, these silicates will vaporize and form an atmosphere whose 
equilibrium vapor pressure can be described by 



surtace 



(7) 



where m is the mass of either a pyroxene or olivine molecule, k is 
Boltzmann's constant , L^ub is the latent hea t of sublimation, and b is 
a constant. We follow iKimura et al] j2002h in our choice of param- 
eters for equation Q. For olivine, m = 169 mn (Mgi iFeo.c)Si04) 
where mn is the mass of atomic hydrogen, L^ub = 3.21 x 10'" erg 
g-', and e* = 6.72 x 10" dyne cm these values were found ex- 
perimentally by evaporating forsterite (Mg2Si04), an end-member 
of olivine (Nagahara et al. 1994^. For pyroxene, m = 60 mn, 
L,ub = 9.61xl0"'ergg~' andg* = 3 . 13x10 " dyne c m ^ these val- 
ues were determined bv lHashimot"ol ( ll990l) for SiO^. lKimura et al] 
( I2OO2I) argue that the parameters for SiOa are appropriate for pyrox- 
ene because evaporation experiments of enstatite (MgSiOs) show 
that SiOa escapes preferenti ally, leaving the miner al with a poly- 
crystalline forsterite coating jTachibana et al.l2002b . We plot equa- 
tion (O for both olivine and pyroxene in Figure [T] 
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According to ISchaefer & Fegle vl ( l2009l) . the atmosphere of a 
hot rocky planet at temperatures similar to Tsurfacc is primarily com- 
posed of Mg, SiO, O, and O2. In our model we choose a mean 
molecular weight for the atmospheric gas of 



// = 30 mn. 



(8) 



similar to the average molecular mass of Mg and SiO. Additionally, 
we set the adiabatic index of the gas to be 7 = cp/cy = 1.3, ap- 
propri ate for diatomic gases at high temperatures (e.g., y^=2000K _ 
1.318. lLange & Forkeill967h . The gas density at the surface of the 
planet is 



(9) 



2.2 Equations Solved 

In our ID hydrodynamic model, we solve the equations for mass, 
momentum, and energy conservation of the gas, assuming a steady 
flow. From mass conservation, 



d , 

^('■Vv) = 0, 
or 



(10) 



where r is the radial distance from the center of the planet, and p 
and V are the density and velocity of the gas. In the frame rotating 
at the orbital angular frequency of the planet, momentum conser- 
vation — with the Coriolis force omitted — implies 



pv 



dv 



dP 

"a7 



GMp _^ 3GMt,pr 



(11) 



where G is the gravitational constant, P = pkT/p is the gas pres- 
sure, T is the gas temperature, and M and M, = OJMq are the 
masses of the planet and the host star (KIC 12557548), respectively. 
The last term on the right-hand side of eq. ( Ill) , which we refer to 
as the tidal gravity term, is the sum of the centrifugal force and the 
gravitational attraction from the star. In deriving eq. (II It , we have 
neglected terms of order {rjaf-. We do not include the contribution 
of the Coriolis force because its magnitude will only be compara- 
ble to the gravitational attraction of the planet when the gas moves 
at bulk speeds comparable to the escape velocity of the planet at 
its Hill sphere (a.k.a. the Roche lobe). This generally occurs near 
the outer boundary of our calculation, which is set by the location 
of the sonic point. Although the sonic point radius and Hill sphere 
radius are distinct quantities, we will find in practice for our mod- 
els that they are close to one another, which is not surprising since 
the wind most easily accelerates to supersonic velocities where the 
effective gravity is weakest — i.e., near the Hill sphere. 

The steady-state energy conservation equation, V • (puv) = 
-P\ ■ V -I- r, together with equation dlOt and the specific internal 
energy u = kT/[(y - l)p], yields 



d 



kT 



(r- l)y" 



kTv dp 

TT + Tcol + Mat- 

u or 



(12) 



The left hand side of equation M2\ tracks changes in internal en- 
ergy, while the terms on the right hand side track cooling due to 
PdV work, heating from gas-grain collisions, and latent heat re- 
leased from grain condensation. We have omitted conduction be- 
cause we found it to be negligible compared to all other terms. 

For gas heating by dust-gas collisions, we assume that each 



collision between a gas molecule and a dust grain transfers an en- 
ergy k{Ti^st - T) to the gas molecule. Then the rate of gas heating 
per unit volume is 



Tcol = -; (rdust - T)p c„ 



(13) 



4.?pi„, p 

where xaustir) = Pdust/p is the spatially varying dust-to-gas mass 
ratio, s = I pm is the assumed grain size, pim = 3 g cm~' is the 
grain bulk density, and 



jkT 



(14) 



is the sound speed of the gas. As discussed in 94.21 our assumption 
that micron-sized grains are entrained in the flow is valid for the 
lowest mass planets we consider (a group that includes, as best we 
will gauge, KIC 1255b) but not the highest mass planets. 
We adopt 



(15) 

a 

for simplicity, so that T^a^t equals T^^^.^^^ at the surface of the planet. 
The dust optical depth to starlight from r to the star is 

3 



4spi, 



■f 



pXi^stdr + T,. 



(16) 



The constant Ts accounts for the optical depth from the outer bound- 
ary of our calculation — the sonic point radius — to the star. We 
arbitrarily set 



T, = 0.01. 



(17) 



Our results are insensitive to our assumed value of t, as long as it 
is « 1. The optical depth at the planet's surface is Tsurf^cc = T(r = 
R), where R is the planetary radius, calculated for a given M by 
assuming a constant bulk density of 5.4 g/cm\ the mean density of 
Mercury. 

Latent heat from grain formation is released at a rate given by 



Tlal - i-subP^'- 



dxd. 



dr 



(18) 



where we have used Ls„b appropriate for pyroxene, as it condenses 
more readily than olivine (see FigurefT] in any case, Ls„(, for olivine 
is not that different). 

Because we do not model dust condensation from first princi- 
ples, the dust-to-gas mass ratio X(ius((r) is a free function, which we 
parametrize as follows: 



l0g[Xdust('-)] 



l0g(j:di,st.max 

1 - arctan 



) - l0g(Xdust,amp) 

S(r-R)'' 



3(r, - R) 



arctan | - 



(19) 



The function Xdust increases monotonically with r to a max- 
imum value of Aidust.max, Starting from a minimum value of 
Jfdi,st.max/i;dusi.amp- We havc two rcasous for choosing Xd„st to in- 
crease rather than decrease with increasing height. The first is phys- 
ical: as the wind launches into space, it expands and cools, allowing 
the metal-rich vapor to more easily saturate and condense (see also 
94.8t . The second is motivated by observations: on length scales 
on the order of the stellar radius, far from the planet, the wind must 
have a relatively high dust content so as to produce the deep transits 
observed by Kepler. 
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As dust condenses from gas and xj„st rises, the gas density 
p must fall by mass conservation. For simplicity we neglect the 
dependence of p on xju,, and so restrict our analysis to xju,, „,ax *k 
1; in practice, the maximum value of x^us, n,ax that we consider is 
10-"^ ~ 0.32. 



Asp, 



(23) 



Note that differences are now computed downwind (Ajx = x,+ i - 
Xj), as V and t have boundary conditions defined at the sonic point 
(j = N). When evaluating variables in (I22t and ( I23t , we average 
across adjacent gridpoints downwind, e.g., p = (pj + pj+i)/2. 



2.3 Relaxation Code 

The structure of the wind is found by the simultaneous solution 
of equations l IlOK 1 11 It . l ll2t . & 1 II6I 1 with the appropriate bound- 
ary conditions. This system constitutes a two-point boundary value 
problem, as the boundary conditions for p and T are defined at the 
base of the flow, while those for v and r are defined at the sonic 
point. 

Although generally more complicated, relaxation codes are 
preferred over shooting methods when solving for the transonic 
Parker wind. This is because for each transonic solution there are 
an infinite number of "breeze solutions" where the bulk speed of 
the wind never reaches supersonic velocities. It is easier to begin 
with an approximate solution that already crosses the sonic point 
and to then refine this solution, than to exhaustively shoot in mul- 
tidimensional space for the sonic point. iMurrav-Clav et all | |2009|) 
used a relaxation code to find the transonic wind from a hot Jupiter. 
Here we follow their methodology; we develop a relaxation code 
based on Section 17.3 of iPress et al., Cl992i) . 



2.3.1 Finite Difference Equations 

We set up a grid of A' = 100 logarithmically spaced radii that run 
from the base of the flow located at the surface of the planet (r = R) 
to the sonic point (r = rj. We replace our system of differential 
equations by a set of firrite-diiference relations. For p and T: 



dp 

/2 ldv\ 



(r- 1) 



T dp LL 

--r + T^CFcoi + r,,,) 

p dr kpv 



(20) 



Ajr = 0, (21) 



where finite differences are calculated upwind (Ajx = Xj - Xj_i), 
as p and T (and by extension £"[ and E2) have boundaiy con- 
ditions defined at the surface of the planet (j = 0). Following 
iPress et al. I( ll992h . we average across adjacent grid points when 
evaluating variables in MQ\ and MVS , e.g., p = (pj-i + P;)/2. The 
finite-difference relations for v and t are: 



E3.J ^ 



V- 



2ykT 7 - 1 



fir 



pv 



(Fcoi + r,a,) 



-GMjr- + 3GM,r/a' 



Ajr = 0, 



(22) 



dr 

E4J = AjT-—Ajr 
dr 



2.3.2 Boundary Conditions 

We need four boundary conditions to solve the system of finite dif- 
ference equations (I20ll-(l23t. At the base of the atmosphere, we set 
the boundary conditions for gas density and temperature using re- 
lations from 92.11 



£^1.0 
£^2.0 



PO - Pvap 
To - T'sui 



■ = 

.ce = 0. 



(24) 



At our outer boundary — the sonic point — we require the bulk 
velocity of gas to equal the local sound speed, and t to equal our 
specified value 02.2t : 



Ei,N 
E4,N 




0. 



(25) 



At every step of the iteration we determine r, by demanding 
dv/dr to be finite at the sonic point. As the denominator of equation 
( 122b vanishes at v = Cj, we require that 



2jkT 7 - 1 



pv 



pr 



(Fcol + Flat) ■ 



GM 3GM,r 



: 0. 



(26) 



2.3.3 Method of Solution 

The relaxation method determines the solution by starting with an 
appropriate guess and iteratively improving it. We use the multi- 
dimensional Newton's method as our iteration scheme, which re- 
quires us to evaluate partial derivatives of all £, j with respect to all 
4A' dependent variables (pj, Tj, Vj, Tj). We evaluate partial deriva- 
tives numerically, by introducing changes of order 10"** in depen- 
dent variables and computing the appropriate finite difi"erences. 
Newton's method produces a 4A' x 4A' matrix, which we invert us- 
ing the numpy library in python. At each iteration we numerically 
solve equation \26\ for r^ and re-map all N gridpoints between R 
and r^. We iterate until variables change by less than one part in 
lO'o. 

We found our iteration scheme to be rather fragile. The code 
only converges if initial guesses are already close to the solution. 
For this reason, we started with a simplified version of the prob- 
lem with a known analytic solution as the initial guess. We then 
gradually added the missing physics to the code, obtaining a solu- 
tion with each new physics input until the full problem was solved. 
Specifical ly, we b egan by solving the isothermal Parker wind prob- 
lem (e.g.. lLamers & Cassinelli 1999). We enforced the isothermal 
condition by declaring 7 to be unity, so that our energy equation 
read dT/dr = (see eq. I2U . Additionally, we set M, and Xdusi to 
zero to remove the efi'ects of tidal gravity and dust. We then slowly 
increased each of the parameters M,, 7, and x^us, (in that order) to 
their nominal values. Any subsequent parameter change (e.g., M) 
was also performed in small increments. 
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Figure 2. Evaporative mass loss rates M vs. planet mass M for isothermal 
dust-free winds. Solid lines are for models that include stellar tidal grav- 
ity, while dashed curves are for models that do not. As the planet mass is 
reduced, all solutions converge toward the free-streaming limit where M 
is not influenced by gravity but instead scales with the surface area of the 
planet (eq. I28t . As inferred from observations, possible present-day mass 
loss rates Mi255b for KIC 1255b are marked in gray. Technically we have 
only a lower limit on Mi255b of O.lMgj / Gyr; for purposes of discussion 
throughout this paper, we adopt 0. l-lMgj / Gyr as our fiducial range (see 
discussion surrounding equation |5). Clearly KIC 1255b cannot be a pure 
pyroxene planet. Subsequent figures will refer to planets with pure olivine 
surfaces (except in j|4.7l where we consider iron). 

In summary, our code contains three main input parameters: 
JCdust.amp and Xdiist.max whlch prescribe the dust-to-gas profile Xi^^i(r), 
and the planet mass M. Our goal is to explore the dependence of 
the mass loss rate 

M = Q.pvr^ (27) 

on these three parameters. Here fl is the solid angle over which 
the wind is launched, measured from the center of the planet; we 
set n = 1 since the high surface temperatures required to produce 
a wind are likely to be reached only near the substellar point. We 
independently vary the parameters Xdu.si,max and Xdu.si,amp to find the 
maximum M for a given M. 



3 RESULTS 

We begin in 33.11 with an isothermal gas model. The isothermal 
model serves both as a limiting case and as a starting point for de- 
veloping the full solution which includes a realistic treatment of 
the energy equation. We provide results for our full model in 33.21 
In the full solution we focus on three possible planet masses M = 
0.01, 0.03, and 0.07 M®, finding that M ^ 0.01 M© yields a maxi- 
mum M that is compatible with the observationally inferred Mi255b. 
We conclude that in the context of our energy equation, the present- 
day mass of KIC 1255b is at most ~0.02 M© (since smaller mass 



planets can generate still higher M that are also compatible with 
^i255b)- In 93.3l we integrate back in time to compute the maximum 
formation mass of KIC 1255b. 



3.1 Isothermal Solution 

We begin by solving a steady, dust-free, isothermal wind described 
by equations l llOt and i ll It with M, set equal to zero. Although this 
is a highly idealized problem, its solution provides us with a starting 
point (i.e., an initial model) from which we are able to solve more 
complicated problems (see g2.3.3t . Our code accepts the standard 
isothe rmal wind solution (e.g.. Chapter 3 of iLamers & Cassinellil 
1 19991) as an initial guess, with minimal relaxation. 

Mass loss rates derived for the isothermal model (7'=2145 K, 
the temperature given by eq. |6]with Tsurface = = 0.01) are plot- 
ted in Figure |2] They depend strongly on planet mass, with large 
gains in M for comparatively small reductions in M, for the basic 
reason that atmospheric densities are exponentially sensitive to sur- 
face gravity. Tidal gravity boosts M by reducing the total effective 
gravity. As M decreases, gravity becomes increasingly irrelevant, 
the thermal speed of the gas eventually exceeds the surface escape 
velocity, and the dependence of M on M weakens. In the "free- 
streaming" limit (dotted lines in Figure |2j, mass loss is no longer 
influenced by gravity, and occurs at a rate 

Mftee ~ PvaporCs^el (28) 

Note that Mfrec decreases when M is reduced, as less surface area 
(R- oc M"'') is available for evaporation. For M ~ O.OSMq, a 
pure olivine surface can produce a (dust-free, isothermal) wind 
for which M becomes compatible with the observationally infeiTed 
value for KIC 1255b, on the order of Mi255b ~ 1 M^/Gyr (this is 
the observed mass loss rate for both dus t and gas combined; se e our 
discussion surrounding equation Island iRappaport et al.ll2012l) . By 
contrast, for a pure pyroxene surface, M «: Mi255b always. Hence- 
forth we will calculate the density Pvapor of gas at the surface using 
parameters appropriate for olivine. 

3.2 Full Solution 

We relax the isothermal condition in our code by slowly increasing 
7 from 1 to 1.3. With this parameter change, we are solving the full 
energy equation (eq.ll2t. Grains are gradually added by modifying 
Xdust.amp and -Xdust.max. necessitating the calculation of optical depth 

In Figure |3] we show the solution for which M was maximized 
(over the space of possible values of Xdust) for M = 0.03 M®, or 
about half the mass of Mercury. Recall from 33.11 that when the 
wind was assumed to be isothermal, a planet of this mass was able 
to reproduce the inferred mass loss rate of KIC 1255b. Now, with 
our treatment of the full energy equation and the inclusion of dust, 
max M plummets by a factor of 40. Relative to the isothermal so- 
lution, the reduction in mass loss rate in the full model is mainly 
caused by the gas temperature dropping from 2095 K at the planet 
surface to below 500 K at the sonic point, and the consequent reduc- 
tion in gas pressure. The temperature drops because gas expands in 
the wind and does PdV work. Gas heating by dust-gas collisions 
or grain formation could, in principle, offset some of the temper- 
ature reduction, but having too many grains also obscures the sur- 
face from incident starlight, decreasing Tsmface and Pvapor (the latter 
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Figure 3. Radial dependence of wind properties in tlie full solution for a 
planet of mass M = O.OSM®. Parameters of the .rjust function are chosen 
to maximize M. In the upper six panels, the full solution is shown by solid 
curves, while the dust-free isothermal solution is marked by dashed curves. 
Values for v have been normalized to the sound speed at 7 = 2145 K, and 
values for p have been normalized to the density po.2i45K = I^Pv^por/ikT) 
evaluated also at T = 2145 K. The x-symbol marks the location of the 
sonic point (the outer boundary of our solution), which occurs close to the 
Hill radius Ruill (marked by a vertical line). The two lower panels show the 
contributions of the individual terms in the energy and momentum equa- 
tions (left and light panels, respectively). 

quantity depends exponentially on the former). The particular dust- 
to-gas profile .Ydust that is used in Figure |3]is such that the ability of 
dust to heat gas is balanced against the attenuation of stellar flux by 
dust, so that M is maximized (for this M). 

Note that the flow at the base begins with velocities of ~10"^^ 
the sound speed. At these subsonic velocities, the atmosphere is 
practically in hydrostatic equilibrium, with gas pressure and gravity 
nearly balancing. Only when velocities are nearly sonic does advec- 
tion become a significant term in the momentum equation (bottom 
right panel of Figure |3). Flow speeds are still high enough at the 
base of the atmosphere to lift dust grains against gravity for all but 
the highest planet masses considered ( 34.21 gas speeds as low as ~1 
m/s in a ~1 yubar atmosphere can enable micron-sized grains to es- 
cape sub-Mercury-sized planets, as can be verified by eauation[30l 
furthermore, the condition that grains not slip relative to gas is con- 
veniently independent of gas velocity in the Epstein free-molecular 
drag regime, as can be seen in equation|29]) 

In Figure|4l we show how M and Tsmfacc vary with the function 



^dust- As just discussed, maxM is reached for a specific choice of 
parameters which balance gas heating by dust and surface obscu- 
ration by dust. (We emphasize, here and elsewhere, that we have 
not identified a physical reason why actual systems like KIC 1255b 
should have mass loss rates M that equal their theoretically allowed 
maximum values.) Apart from the region near the peak in M, con- 
tours of constant Tsurfjicc are roughly parallel to contours of constant 
M, suggesting that the value of Tsmfacc is more important for deter- 
mining M than the specific functional form of jc^us,. This finding 
increases the confidence we have in the robustness of our solution. 

It is clear from Figures |3]and|4]that a planet of mass M = 0.03 
M0 can only emit winds for which max M < min Mi255b ~ 0. IM0 / 
Gyr (at least within the context of our energy equation). In Figure|5] 
we show results for M = 0.01 M^, for which max M > min Mnssb- 
For this M, the mass loss rate M is maximized when the planet 
surface is essentially unobscured by dust. At our arbitrarily chosen 
value for Xd„st ~ 3x 10"^, the wind is not significantly heated by dust 
and expands practically adiabatically. The maximum M is reached 
for an essentially dust-free solution because the wind is blowing at 
too high a speed for dust-gas collisions to be important. That is, the 
timescale over which gas travels from the planet's surface to the 
sonic point is shorter than the time it takes a gas particle to collide 
with a dust grain: gas is thermally decoupled from dust. Were we to 
increase xjas, to make dust-gas heating important, the flow would 
become optically thick and the wind would shut down. If the model 
shown in Figure|5]does represent KIC 1255b, then the dust grains 
that occult the star must condense outside the sonic point, beyond 
the Hill sphere. 

Figure |6] shows the solution for M = 0.07 Mgj. In this case M 
is maximized when enough dust is present to heat the gas signifi- 
cantly above the adiabat. At this comparatively large planet mass, 
initial wind speeds are low enough that a gas particle collides many 
times with dust over the gas travel time, so that dust-gas collisions 
heat the gas effectively near the surface of the planet. Further down- 
wind, near the sonic point, latent heating overtakes PdV cooling 
and actually increases the temperature of the gas with increasing 
altitude. The wind is much more sensitive to heating near the base 
of the flow than near the sonic point: if we omit latent heating from 
our model — so that gas cools to ~500 K at the sonic point — then 
M is reduced by only a factor of two. In comparison, M drops by 
several orders of magnitude if dust-gas collisions are omitted. 

The dependence of M on x^^^t in both the "dust-free low- 
mass" and the "dusty high-mass" limits is illustrated in Figure |7] 
In the low-mass limit, M = max M when no dust is present. In this 
regime, gas moves too quickly for heating by dust to be significant; 
dust influences the flow only by attenuating starlight, and as long 
as Tsurfacc <K 1, dust hardly affects M. By contrast, in the high-mass 
limit, M = max M for the dustiest flows we consider (i.e., the high- 
est values of Jfdust.max)- In this regime, M is very sensitive to dust 
abundance, with values spanning six orders of magnitude over the 
explored parameter space. 

In Figure [8] we show mass loss rates as a function of planet 
mass for the full model. We emphasize that these are maximum 
mass loss rates, found by varying Xd„st- At low planet masses, 
maxM for the full model converges with M for the isothermal 
model because both approach the free-streaming limit, where grav- 
ity becomes irrelevant and M is set entirely by conditions at the 
surface of the planet (eg. 128). To reach Mi255b > 0. IMgj / Gyr, the 
present-day mass of KIC 1255b must be < 0.02 Mgj, or less than 
about twice a lunar mass. 
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Figure 4. Left panel: Dependence of M on dust abundance for a planet of mass M = O.OSMq. The dust abundance parameters JCdust, miix ™d JCd^st, amp ar£ 
defined in equation U9t : see also Figure|3] Both parameters are varied independently to find the maximum M, which occurs at log(X;jus( j„„p, ;Cdust,max) ~ (1-7. 
-2.2). The full model con'esponding to this maximum M is detailed in Figure|3] Each black dot corresponds to a full model solution. Colour contours for M 
are interpolated using a cubic polynomial. Right panel: Optical depth Tsurfa^c between the star and the planetary surface as a function of dust abundance for a 
planet of mass M = O.OSMgj. Apart from the region near the peak in M, contours of constant Tsurhce are roughly parallel to contours of constant M, suggesting 
that the value of Ts^pface is more important for determining M than the specific functional forni of Ji,iust. 



We have verified a posteriori for the models shown in Figures 
[3] [5] and|6]that the sonic point is attained at an altitude where the 
colhsional mean free path of gas molecules is smaller than r^, so 
that the hydrodynamic approximation embodied in equations lllOt- 
l ll2t is valid. In other words, the exobase lies outside the sonic point 
in these M = max M models. The margin of safety is largest for the 
lowest mass models. 



3.3 Mass-Loss Histories 

We calculate mass-loss histories M{t) by time-integrating our solu- 
tion, shown in Figure [8] for max M{M). Before we integrate, how- 
ever, we introduce < f^^y < 1 to account for the duty cycle of 
the wind: we define /d„ty ■ max M as the time-averaged mass-loss 
rate. We estimate that /d^y ~ .5 bas ed on the statistics of transit 
depths compiled by Brogi et al. l l2012l see their figure 2). We time- 
integrate /duty • max M to obtain the Mil) curves shown in Figure|9] 
(This factor of 2 correction for the duty cycle should not obscure 
the fact that the actual mass loss rate, time-averaged or otherwise, 
could still be much lower than the theoretically allowed max M, a 
possibility we return to thoughout this paper.) 

We highlight the case of a planet with a lifetime of f^fe = 5 
Gyr. Under our full (non-isothermal) model (right panel of Figure 
|9]l, such a planet has a mass at formation of 0.06 Mgj and grad- 
ually erodes over several Gyr until it reaches M ~ 0.03 — 
whereupon the remaining mass is lost catastrophically over a short 
time (just how short is estimated for the specific case of KIC 1255b 
below). Contrast this example with planets having initial masses 



> 0.07 M0 — these have such low M that they barely lose any 
mass over Gyr timescales. 

For comparison, we also present mass-loss histories using a 
dust-free isothermal model with stellar tidal gravity (left panel of 
Figure |9). Compared to our non-isothermal solutions — for which 
gas temperatures fall immediately as the wind lifts off the planet 
surface — the isothermal solution corresponds to a flow which 
stays relatively pressurized and which therefore enjoys the largest 
M for a given M. The isothermal wind thus represents an endmem- 
ber case. However, the mass-loss histories under the isothermal ap- 
proximation do not differ qualitatively from those using the full 
energy equation. For isothermal winds, the dividing mass between 
planet survival and destruction within 10 Gyr is about 0.11 M®, 
only ~40% larger than the value cited above for our non-isothermal 
solutions. 



3.4 Possible Mass Loss Histories for KIC 1255b 

Figures [8] and |9] can be used to sketch possible mass-loss histories 
for KIC 1255b. For our first scenario, we employ the full model 
that solves the full energy equation. Today, to satisfy the obser- 
vational constraint that Mi255b > 0. 1 M0/Gyr (see the derivation 
of our equation |5}, the planet must have a present-day mass of at 
most M X 0.02 Mgj, or roughly twice the mass of the Moon (Fig- 
ure [8}. Such a low mass implies that currently KIC 1255b is in a 
"catastrophic evaporation" phase (vertical straight lines in Figure 
O. Depending on the age of the planet (i.e., the age of the star: 
1-10 Gyr), KIC 1255b originally had a mass of at most 0.04-0.07 
M®: 2-4 times its maximum current mass, or about the mass of 
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Figure 5. Same as Figure|3] but for a planet of mass M = O.OIM^. For this 
planet mass, M is maximized when essentially no dust is present (i.e., the 
atmosphere is essentially transparent) as gas moves too quickly for dust-gas 
collisions to heat the gas. As such, all heating tenns are negligible and the 
gas expands adiabatically. 
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Figure 6. Same as Figure[3] but for a planet of mass M = O.OTMg. For this 
mass, M is maximized for the dustiest flows considered in this paper The 
outflow is launched at such low velocities that coUisional dust-gas heating 
is important. Near the sonic point, latent heating overtakes PdV coohng and 
the temperature of the gas rises. 



Mercury. Starting from today, the time the planet has before it dis- 
integrates completely is on the order of fevap ~ M/ (/duty • max M), 
which is as long as 400 (40) Myr for M = 0.02 M®, f^^ty = 0.5, 
and max M = 0.1 (1) M® / Gyr. 

If our energy equation were somehow in en'or and the wind 
better described as isothermal, then the numbers cited above would 
change somewhat. The present-day mass of KIC 1255b would be 
at most M ~ 0.01 Mq; the maximum formation mass would be be- 
tween 0.08-0.11 M®; and the evaporation timescale starting from 
today could be as long as 4vap ~ 400 Myr if the present-day mass 
M ~ 0.07M©. 

As we have stressed throughout, a critical assumption we have 
made is that the time-averaged mass loss rate is /duty • max M with 
/duty as high as 0.5. We have not identified a physical reason why 
the actual mass loss rate should be comparable to the theoretically 
allowed max M for a given M. If our assumption were in error and 
actual mass loss rates M <c max M, then the present-day mass of 
KIC 1255b, the corresponding evaporation time, and the formation 
mass of KIC 1255 would all decrease from the upper bounds cited 
above. 

The degeneracy of possible masses and evaporative histories 
outlined in this section could be broken with observational searches 
for the progenitors of KIC 1255b-like objects ( 94.61 . 



4 DISCUSSION 

We discuss how mass loss is not "energy-limited" in 94.11 how 
gas entrains dust in 94.21 how the evaporative wind can be time- 
variable in 94.31 the effects of winds on orbital evolution in 94.41 
the interaction of the planetary wind with the stellar wind in 94.51 
the occurrence rate of quiescent progenitors of catastrophically 
evaporating planets in 94.61 the possibility that evaporating planets 
can be naked iron cores in 94.71 and the ability of dust to condense 
out of the wind in 94.81 



4.1 Mass Loss is Not Energy-Limited 

Atmospheric mass loss rates are commonly assumed to be "UV- 
energy-limited": it is assumed that a fixed, order-unity fraction e of 
the incident UV flux Fuv does PdV work and lifts material out of 
the planet's gravitational well. Under the energy-limited assump- 
tion, the mass loss rate is M ~ enR'Fuy /(GM/R). T he validity of 
this formula has be en tested by Watson et al. I( ll98lh for terrestrial 
atmospheres and bv lMurrav-Clav et al.l ( l200S ) for hot Jupiters. 

The energy-limited formula does not apply at all to the evap- 
orating rocky planets considered in this paper. The stellar UV flux 
is not essential for low-mass planets because their escape veloci- 



© 0000 RAS, MNRAS 000, 000-000 



10 D. Perez-Becker (ir E. Chiang 




Figure 7. Left panel: Same as the left panel of Figure|4] but for a planet of mass M = 0.006Mg). In this low-mass limit, M is maximized when no dust is 
present (i.e., when the atmosphere is transparent) because gas is moving too quickly for heating by dust to be significant. Right panel: Same as left panel, but 
for a planet of mass M = O.lMgj. In this high-mass limit, M is maximized for the dustiest flows we consider (i.e., the highest values of .Vdust.max) because flow 
speeds near the surface are slow enough for dust-gas energy exchange to be significant. 



ties are so small that energy deposition by photoionization is not 
necessary for driving a wind. Even optical photons — essentially 
the stellar bolometric spectrum — can vaporize close-in planets. 
Merely replacing Fuv with Fboiomctric in the energy-limited formula 
would still be misleading, however, because most of the incident 
energy is radiated away and does no mechanical work. More to the 
point, e is not constant; as discussed in 33.21 the energetics of the 
flow changes qualitatively with planet mass. If we were to insist 
on using the energy-limited formula, we would find e ~ 10"* for 
M ~ O.IM®, and e ~ 10"^ for M ~ O.OIM© (see Figurei}. 



4.2 Dust-Gas Dynamics 

In our model we have assumed that dust grains condensing within 
the wind are carried along without any relative motion between dust 
and gas (1 -fluid approximation). Grains are well entrained if their 
momentum stopping times 

'"grain Vrel .Spi„, 

'stop (^y) 

are shorter than the grain advection time fadv ~ ^/i'- In equation 
([29}, 

'^^giain is the mass of an individual grain and v,-ei is the relative 
gas-grain velocity. The aerodynamic drag force in the denominator 
is given by the Epstein law, which is appropriate for I'lei < Cs and 
grain sizes s smaller than the gas coUisional mean free path (for our 
flows, /Imfp ~ fi/(pcr) > 2 cm, where cr ~ 3 x 10"'^ cm^). 

We find using our full model (for M = max M) that fstop/fadv < 
1 everywhere for M < 0.03 M®, confirming our assumption that 
micron-sized (and smaller) grains are well-entrained in the winds 
emanating from such low-mass planets. When M = 0.03 Mq, 
fstop < 4dv close to the planet's surface, but as gas nears the sonic 



point, fstop becomes comparable to 4dv For M > 0.03 M0, the 1- 
fluid approximation breaks down near the sonic point. For such 
large planet masses, future models should account for gas-grain 
relative motion — in addition to other effects that become increas- 
ingly important near the sonic point / Hill sphere boundary (e.g., 
Coriolis forces and stellar radiation pressure). Note that the 1-fluid 
approximation should be valid for the present-day dynamics of KIC 
1255b, since according to the full model its cuiTent mass is at most 
M « 0.02 M® < 0.03 M®. 

Although dust can slip relative to gas, dust can still be trans- 
ported outward and escape the planet if the aerodynamic drag force 
exceeds the force of gravity: 

PVrelCsi^ >Pmt^^g, (30) 

where g = G(M/r^ - 3M,r/a') is the total gravitational accelera- 
tion. Grains that manage to be lifted beyond the Hill sphere (which 
is situated close to the sonic point in all our models) are no longer 
bound to the planet. Outside the Hill sphere, dust decouples from 
gas and is swept into a comet-like tail by stellar radiation pressure 
and Coriolis forces. 

For M « 0.03 M® and our adopted grain size .s = 1 yum, 
the forces in equation ( I30t are comparable near the sonic point 
when I'l-ei ~ Cs (i.e., when grains are barely lifted by drag). For 
M > 0.03 M®, micron-sized and larger grains will not be dragged 
past the sonic point according to i30\ — thus our assumption that 
they do fails for such high-mass planets, and will need to be recti- 
fied in future models. Smaller, sub-micron-sized grains can, how- 
ever, be lifted outward. Moreover, sub-micron sized grains may 
heat the gas more effectively because they have greater geomet- 
ric surface area for dust-gas collisions (at fixed Xdu.si) and because 
their temperatures exceed those of blackbodies (their efficiencies 
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Figure 8. Maximum mass loss rates M vs. planet mass M for the full model. 
At each x-marked mass, max M was found by varying .vj„st as described, 
e.g., in Figure |4] The solid curee is a cubic spline interpolation. Mass loss 
rates for the full model are generally lower than for the isothermal model 
(dashed curve) because the dust-gas heating terms we have included in our 
full model turn out to be inefficient. At low planet masses, the isothermal 
and full models converge because both approach the free-streaming limit, 
where gravity becomes irrelevant and M is set entirely by conditions at 
the surface (i.e., surface area, equilibrium vapor pressure, and sound speed; 
eg. 1281 . According to the full model, the present-day mass of KIC 1255b 
is < 0.02 Mgj, or less than twice the mass of the Moon; for such masses, 
max M > min Mnssb ~ 0. 1 / Gyr 
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Figure 9. Mass-loss histories M(f) obtained by time-integrating /jmy ■ 
max M(M) with /j^ty = 0.5 for the isothermal and full models. We high- 
light the case of a planet with a 5 Gyr hfetime. For the full model (right 
panel), this corresponds to an initial mass of ~0.06 M^, slightly larger than 
the mass of Mercuiy. Such a planet could have formed in situ and slowly 
eroded over several Gyr until reaching ~0.03 Mg;, whereupon the planet 
evaporates completely in a few hundred Myr By contrast, planets with for- 
mation masses ^ O.OVMgj survive (in the full model) for tens of Gyrs with- 
out significant mass loss. 



for emission at infrared wavelengths are much less than their ab- 
sorption efficiencies at optical wavelengths). The superior momen- 
tum and thermal coupling enjoyed by grains having sizes .s <k 1 jim, 
plus their relative transparency at optical wavelengths — which 
helps them avoid shadowing the planet surface from stellar radi- 
ation — motivate their inclusion in the next generation of models. 

Our model could be improved still further by self-consistently 
accounting for how the gas density p must decrease as dust grains 
condense out of gas. Moreover, the restriction that Xdust < 1 could 
be lifted. 



4.3 Time Variability 

Occultations of KIC 12557548 vary in depth from a maxi- 
mum of 1.3% to a m inimum of < 0.2% on orbital timescales. 
iRappaport et al.l ( l2012h discussed qualitatively the possibility that 
such transit depth variations arise from a limit cycle that alter- 
nates between high-M and low-M phases. A high-M phase that 
produces a deep eclipse would also shadow the planet surface from 
starlight. The resultant cooling would lower the surface vapor pres- 
sure and lead to a low-M phase — after which the atmosphere 
would clear, the surface would re-heat, and the cycle would begin 
anew. Such limit cycle behavior could be punctuated by random 
explosive event s that release dust, similar to thos e observed on lo 
I Geisslej|2003l ; see also the references in §4.2 of IRappaport et al] 
l2012h . 

Order-of-magnitude variations in M arise from only small 
fractional changes in surface temperature because the vapor pres- 
sure of gas over rock depends exponentially on temperature. For 
example, from equations (|6) and (|7), we see that increasing the sur- 
face optical depth Tsuriace from 0. 1 to 0.4 reduces the surface temper- 
ature Tsui-facc by 150 K and the base density pvapoi by a factor of 10. 
Such changes would cause M to drop by more than a factor of 10, 
because M scales super- linearly with Pvapor: a linear dependence re- 
sults simply because M scales with gas density, while an additional 
dependence arises because the wind speed increases with the gas 



To reproduce the orbit-to-orbit variations in transit depth ob- 
served for KIC 1255b, the dynamical time fdyn of the wind cannot 
be much longer than the planet's orbital period of P^rb = 15.7 hr. 
The dynamical time is that required for dust to be advected from 
the planet surface to the end of the comet-like tail that occults the 
star: it is the minimum timescale over which the planet's transit 
signature "refreshes". If fjyn 3> P^rh, then we would expect transit 
depths to be correlated from one orbit to the next — in violation of 
the observations. 

Referring to the full model for a KIC 1255b-like mass of M x 
O.OIM0 (Figure|5]l, we estimate that 



hr + 



OAR, 



14 hr, 



(31) 



where we have split f^y,, into two parts: the first integral is the time 
for dust to reach the sonic point (the outer boundary of our calcu- 
lation), while the second integral is the time for dust to travel out 
to 0. IR, (a circular, optically disk of this radius would generate a 
1% transit depth). The first integral is performed numerically us- 
ing our full model (upper left panel of Figure |5), while the second 
integral is estimated to order-of-magnitude using a characteristic 
grain velocity (well outside the planet's Hill sphere) of vi ~ 10 km 
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s"' ( IRappaport et al.ll2012l : see their equation 6 and related com- 
mentary). Since ?dyii £ ^'orb, we conclude that the wind/cometary 
tail can refresh itself quickly enough to change its appearance from 
orbit to orbit. 

We can go one step further. The timescale over which M 
changes should be the timescale over which the stellar insolation 
at the planetary surface changes — in other words, the timescale 
over which the "weather" at the substellar point, where the wind 
is launched, changes from, e.g., "overcast" to "clear". Ignoring 
the possibility of volcanic eruptions, we estimate this variability 
timescale as the time for the wind to reach the Hill sphere bound- 
ary, at which point the Coriolis force has turned the wind by an 
order-unity angle away from the substellar ray joining the planet to 
the star (i.e., beyond the Hill sphere, the dust-laden wind no longer 
blocks stellar radiation from hitting the substellar region where the 
wind is launched). Because the Hill radius is situated near the sonic 
point, this variability timescale is given approximately by the first 
integral in equation ( I3U : f^^^ ~ 13 hr. 

The fact that f, is neither much shorter than nor much longer 

dyn o 

than Porb supports our proposal that the observed time variability of 
KIC 1255b is driven by a kind of limit cycle involving stellar inso- 
lation and mass loss. Clearly if /^^^^ » Pa^^,, orbit-to-orbit variations 
in M would be impossible. Conversely, if t^^^ « Porb — or more 
precisely if f^^^_^ were much shorter than the transit duration of 1.5 
hr — then the wind would vary so rapidly that each transit observa- 
tion would time-integrate over many cycles, yielding a smeared-out 
average transit depth that would not vary from orbit to orbit as is 
observed. 



4.4 Long-Term Orbital Evolution 



When computing planet lifetimes in 33.31 we have assumed that the 
planet does not undergo any orbital evolution while losing mass. 
This approximation is valid because gas leaves the planet's surface 
at velocities viaunch ^ Q ~ 1 km s"'. Launch velocities are so low 
compared to the orbital velocity of the planet Vo^b ~ 200 km s"' that 
the total momentum imparted by the wind is a tiny fraction of the 
planet's orbital momentum. 

We quantify this as follows. We estimate the total change in 
sem imajor axis a and eccentr icity e from Gauss' equations (see, 
e.g., lMurrav & DermotlllT999l) . Consider the case in which the gas 
is launched at a velocity via^nch at an angle a from the substellar ray 
and in the plane of the planet's orbit. The angle a could be non- 
zero because of surface inhomogeneities or asynchronous rotation 
of the planet. Neglecting terms of order e' and ae, we find that a 
and e evolve according to 



da 
dt 
de 

Tt 



I'kimch M 

a (e + a) , 

Vorb M 
Vlaunch M 



Vorb M 



(l+a), 



(32) 



where the first term in parentheses on the right-hand-side of each 
equation is due to the radial component of the perturbation and the 
second term is due to the azimuthal component (where radius and 
azimuth are measured in a cylindrical coordinate system centered at 
the star and in the plane of the planet's orbit). Integrating equations 
( I32t yields the changes in a and e accumulated over the planet's 
age: 









1 ^'launch 


a 




I, ^filial 




Ae - 


- M 


Mnlial \ 


^'launch 






Mfinal / 


^'oi-b 



(33) 



When evaluating equation ( I33t . we have used the maximum possi- 
ble launch velocity v'launch ~ Cs (valid in the "free-streaming" limit) 
and M,„,„ai/Mfl„ai ~ 5 (see 

The fact that evaporating planets undergo negligible orbital 
evolution suggests they could have resided on their cu rrent orbits 
for Gy rs — indeed they might even have formed in .vi'rif. lSwift et al.l 
( 120131) would disfavor in-situ formation of KIC 1225b because at 
the planet's orbital distance, dust grains readily sublimate, and os- 
tensibly there would have been no solid material in the primor- 
dial disk out of which rocky planets could have formed. How- 
ever, a less strict in-situ formation scenario is still viable. Solid 
bodies larger than dust grains obviously have longer evaporation 
times. Such large planetesimals could have drifted inward through 
the primordial gas disk and then assembled into the progenitors 
of objects like KIC 1 255b at their cun 'ent close-in distances (e.g^ 
l Youdin & Shm200a : (Hansen & Murray 2012: Chiang & Laughliiil 
l2012h . Solid particles can avoid vaporization if they merge faster 
than they evaporate. 



4.5 Flow Confinement by Stellar Wind 

When computing the transonic solution for the planetary wind, we 
implicitly assumed that the outflow was expanding into a vacuum. 
In reality, the outflow from the planet will collide with the stellar 
wind and form a "bubble" around the planet. For our solution to be 
valid, the radius of the bubble — i.e., the surface of pressure bal- 
ance between the winds — has to be downstream of the sonic point 
where the flow is supersonic. This way, any pressure disturbance 
created at the interface cannot propagate upstream and influence 
our solution. 

At the wind-wind interface, the normal components of the 
pressures of the planetary and stellar winds balance. The total pres- 
sure includes the ram, thermal, and magnetic contributions. We 
compare the pressure of the stellar wind at a ~ O.OI AU ~ 4R. 
with that of the planetary wind computed at to gauge whether 
the interface occurs downstream of r^. For the pressure P, of a 
main-sequence, solar-type star we are guided by the Solar wind. 
Flow speeds of the Solar wi nd are measured by tracking the tra- 
jector ies of coronal features JSheelev et al.lll997t lOuemerais et al.l 
l2007h . At a = O.OI AU the Solar winc||is still accelerating with typ- 
ical flow speeds of v, ~ 100 km s"' . For this v,, and a canonical So- 
lar mass-loss rate of 2 x lO"'"* Mq yr"' , the proton number density is 
n, = 2x10^ cm~^.W e take the proton temperature to be T, ~ 10'' K 
( ISheeleYetalJl997l) and the heliospheric magnetic field strength to 



^ These stellai' wind parameters are appropriate for the "slow" Solai' wind 
blowing primarily near the equatorial plane of the Sun. In addition, the Solar 
wind also contains a "fast" component which emerges from coronal holes. 
During Solar minimum, the fast Solar wind is confined mainly to large 
heliographic latitudes, b ut may reach the equator plane during periods of 
increased solar activity l lKohl et al.|[T998l : iMcComas et"ai]|20O3l) . Because 
most planets are expected to orbit near their stellar equatorial planes, we 
take the slow component of the Solar wind to be a better guide. 
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be B. « 0.1 G at 0.01 AU jKim et alj|2012h . The total stellar pres- 
sure at 0.01 AU is then P, ~ n.mn\i + n.kT. + Bl/(Sn) ~ 4 x lO""* 
dyne cm"^, with the magnetic pressure dominating other terms by 
an order of magnitude. 

The ram and thermal pressures of the planet's wind are equal 
at the sonic point and add up to (for M = maxM) ~ 2 x 10"' 
(1 X 10"'*) dyne cm"^ when M = 0.03 (0.07) M^. Thus for M < 
0.03 Mq (a range that includes possible present-day masses of 
KIC 1255b), Ps > P, so that the planetary wind blows a bubble 
that extends beyond the sonic point and the transonic solutions we 
have computed are self-consistent. By contrast for M ~ 0.07 M^, 
Ps ~ P,/4 < P, so that the stellar wind pressure will balance the 
planetary wind pressure inside of rj. The stellar wind will prevent 
the planetary wind from reaching supersonic velocities; the plan- 
etary outflow will conform to a "breeze" solution with a reduced 
M. If the colliding winds in our problem behave like those of hot 
Jupiters and thei r host stars, then Ps ~ P,/ 4 will reduce M by ~30% 
(see figure 12 of lMurrav-Clav et al.ll2009l) . 

We have found P, to be dominated by magnetic pressure. 
Since magnetic fields only interact with the ionized component of 
the planetary wind, and since the planetary wind may not be fully 
ionized, we might be over-estimating the efi'ect of P,. Dust grains in 
the planetary wind may absorb free charges so that the wind may be 
too weakly ionized to couple with the stellar magnetic field. If we 
ignore the magnetic contribution to the stellar wind, then P^ > P, 
for all planet masses that we have considered, and all our transonic 
solutions would be self-consistent. 

What about the planet's magnetic field? A sufficiently ion- 
ized outflow could be confined by the planet's own magnetic field. 
To assess the plausibility of this scenario, we use Mercury's field 
as a guide for KIC 1255b. The flyby of Mariner 10 and recent 
measurements by MESSENGER show that Mercury possesses a 
weak magnetic field with a surface stre ngth of ~ 3 x 10"' Gauss 
jNess etal.|[T974l :l Anderson et al. I l2008h . The associated magnetic 
pressure ~ 4 x 10"^(/?/r)* dyne cm"^ is small compared to both the 
thermal pressures that we have computed at the planet surface and 
the hydrodynamic pressures at the sonic point. For planetary mag- 
netic fields to confine the wind, they are required to have surface 
strengths in excess of ~ 30 Gauss. 



4.6 Occurrence Rates of Close-in Progenitors 

According to our analysis in 33.41 cuiTently KIC 1255b has a mass 
of at most 0.02-0.07 Mq, and is in a final, short-lived, catastroph- 
ically evaporating phase possibly lasting another /evap ~ 40-400 
Myr during which dust is ejected at a large enough rate to produce 
eclipse depths of order 1%. But for most of its presumably Gyrs- 
long life, KIC 1255b was a more quiescent planet — larger but still 
less than roughly Mercury in size ( g3.4b , with a stronger gravity 
and emitting a much more tenuous wind. We therefore expect that 
for every KIC 1255b-like object discovered, there are many more 
progenitors in the quiescent phase — i.e., planets having sizes up to 
about that of Mercury orbiting main sequence K-type stars at ~0.01 
AU. 

Out of the ~ 45,000 main sequen ce K-type stars in the Kepler 
Input Catalogue teatalha et all 1201 Ol) . there is apparently only 1 
object like KIC 1255b. Then 

yobsei ved ^rogenitoi ./transit /evap . cnr\n ' (3^) 



103 p ' M 




Figure 10. Mass loss rates M, computed using the isothennal model atT = 
2145 K, for iron planets and olivine planets of varying M. In the low-M, 
free-streaming limit, the higher vapor pressure of iron allows for a higher 
M than for olivine. At lai'ger M, mass loss rates are lower for iron planets 
than for olivine planets because of the higher molecular weight of iron. Our 
estimate for the maximum present-day mass of KIC 1255b varies by a factor 
of two between the iron and olivine scenarios. 

where /pioecnitoi is the intrinsic occurrence rate of progenitors around 
K stars with orbital periods shorter than a day, /transit ~ RJa ~ 1/4 
is the geometric probability of transit, and /evap = fevap/fiife ~ 
(40, 4001 Myr / 5 Gyr ~ (0.8,8)% is the fraction of the planet's 
lifetime spent in the catastrophic mass-loss stage. Inverting equa- 
tion l |34| (, we estimate that /p,ogcnitoi ~ ( 1, 0. 1 }% of K stars harbor a 
close-in planet having less than the mass of Mercury. 

For each KIC 1255b-like object there should be /^^[^ = 
fiife/fevap ~ (130, 13) plaucts with radii no larger than about Mer- 
cury's (< OARq), transiting K-type stars with sub-day periods. 
These planets have transit depths of (R/R,)- < 10"^ — small 
but possibly detectable by Kepler if light curves are folded over 
enough periods (S. Rappaport 2012, personal communication) and 
if M ~ max M so that the progenitor masses attain their maximum. 
Mercury-like values. If M were actually maxM, the progen- 
itor sizes would be less than that of Mercury. The smaller sizes 
would decrease their lifetimes /evap and thus increase their expected 
number /„ap. but would also render them undetectable even with 
Kepler. 

4.7 Iron Planets 

In our model, we found the (maximum) present-day mass of KIC 
1255b to be about 1/3 of its (maximum) mass at formation. If KIC 
1255b began its life similar in composition to the rocky planets in 
the Solar System, evaporation may have stripped the planet of its 
silicate mantle, so that only its iron core remains: KIC 1255b could 
be an evaporating iron planet today. 
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We estimate mass loss rates of a pure iron planet using our 
isothermal model, including tidal gravity. We set the mean molec- 
ular weight of the gas to /jpe ~ 56/nH and use a bulk density for 
the plane t of 8.0 g cm~^, appr opriate for a pure iron planet of mass 
O.O1M0 l lFortnev et al.[ 20071) . We fit lab oratory meas urements of 
the iron vapor pressur^ at T = 2200 K jDesai[l98^ to equation 
(|7j, obtaining = 7.8 x 10" dy ne cm~^ for a latent heat of subli- 
mation Lsi,b = 6.3 X 10'" erg g"' toesai|[l983) and m = /jpc- 

In Figure [10] we compare mass loss rates derived for the 
isothermal model with T = 2145 K for both an iron planet and a 
pure olivine planet. On the one hand, at this temperature, the vapor 
pressure for iron Pvapor.Fe = 1-8 x 10^ dyne cm"^ is about 50 times 
higher than for olivine, which raises the mass loss rate by a similar 
factor in the free-streaming limit appropriate for low masses (see 
equation[28]). On the other hand, at higher masses, M for iron drops 
below that for olivine because the iron atmosphere, with its higher 
molecular weight, is harder to blow off. These two effects counter- 
act each other so that an iron planet and a silicate planet both reach 
^i255b ~ 1 ^© / Gyr at a similar M. Thus at least within the con- 
text of isothermal winds, our estimate for the present-day mass of 
KIC 1255b is insensitive to whether the evaporating surface of the 
planet is composed of iron or silicates. Figure [TTIdisplavs mass loss 
histories for an iron planet using our isothermal model at T = 2145 
K. Initial planet masses are up to a factor of two lower than those 
of our full model using silicates. 

According to Figure [TT] iron planets with M > O.O5M0 
survive for tens of Gyrs. Compare this result to its counterpart 
in Figure |9] (left panel), which shows that olivine planets with 
M < 0.1 M0 evaporate within ~10 Gyr. This comparison suggests 
that for planets with the right proportion of silicates in the mantle 
to iron in the core, mantles may be completely vaporized, leaving 
behind essentially non-evaporating iron cores. Such massive, qui- 
escent iron cores might be detectable by Kepler via direct transits 
(see 34.6I 1. 



4.8 Dust Condensation 

Recall that we have not modelled the microphysics of dust forma- 
tion, but treated the dust-to-gas ratio as a free function. Are the 
conditions of our flow actually favorable for the condensation of 
dust? Generally, any gas whose partial pressure is greater than its 
vapor pressure (at that T) will condense and form droplets. Conden- 
sation proceeds until all vapor in excess of saturation is in cloud 
particles. If the condensates have sedimentation velocities larger 
than updraft speeds they will ra in out of the atmosphere; other- 
wise they remain al oft (see. e.g..lLewisl | l969l; iMarlev et al .111 9991 ; 



0.07 



lAckerman & Marlevll200ll ; lHelling et al.ll200lh 

In Figure [T2I we compare gas pressures of our solutions (for 
M = maxM) at M = 0.01, 0.03, and 0.07 M@ with the equilibrium 
vapor pressures of olivine and pyroxene. At the base of the flow, 
the gas pressure equals the saturation vapor pressure of olivine by 
construction. As gas expands and cools, its pressure generally re- 
mains above the saturation pressure of pyroxene, so conditions are 
favorable for the condensation of pyroxene grains. Unfortunately 
we have an embarrassment of riches: at T < 1700 K, our gas pres- 
sures are many orders of magnitude higher than silicate vapor pres- 



' We have found measurements of the iron vapor p ressure at T ~ 20 00 K 
to differ by a factor of two in the literature (see, e.g.. lNuth et al .l2003h . 




Figure 11. Mass-loss histories M{1) for an iron planet, obtained by time- 
integrating /juty ■ max M with /j^ty = 0-5 for an isothermal wind. An iron 
planet with a 5 Gyr lifetime will have an initial mass of ~0.044 M^. As 
was the case for olivine planets (see Figure|9), the catastrophic evaporation 
stage lasts only for ~100 Myr. Iron planets (or planetary iron cores) with 
masses greater than 0.05 survive for over 10 Gyr. 



sures, and the concern is that gas condensation will lower pres- 
sures precipitously to the point where winds shut down. This fate 
could be avoided if grains take too long to condense out of the 
outflowing and rapidly rarefying gas — i.e., the finite timescale 
of grain condensation might permit gas to remain supersaturated. 
Furthermore, even before bulk condensation can begin, seed parti- 
cles must first nucleate from the gas phase. These issues call for a 
time-dependent analysis of grain formation, perhaps along the lines 
made for brown dwarf and exoplanet atm ospheres jllelling et al.l 
l200lLl2008al lbl; lHelling & Rietmeiieil2"009l) . Accounting for the ex- 
tra heating from small grains with sizes <k 1 /jm as they first con- 
dense out of the wind ( 94.21 would keep the gas closer to isothermal 
and help to prevent catastrophic condensation. 

Of course, by the wind leaves the sonic point and is diverted 
into the trailing comet-like tail, it must be full of grains to explain 
the observed occultations of KIC 1255b. In this paper, we have not 
modeled the dusty tail at all. Future work should address the dy- 
namics of the tail and compute its extinction profile, with the aim 
of reproducing the observed transit lightcurves. 



5 CONCLUSIONS 

Our work supports the hypothesis that the observed occultations of 
the K-star KIC 12557548 originate from a dusty wind emitted by 
an evaporating planetary companion ("KIC 1255b"). We reach the 
following conclusions. 

(i) Maximum present-day mass 

Our estimates for the maximum present-day mass of KIC 1255b 
range from max M x Q.QIM^ (roughly twice the mass of the 
Moon) to max M ^ 0.07 (slightly larger than the mass of Mer- 
cury), depending on how strongly the wind is heated by dust and 
can remain isothermal (Figure|8). For these maximum masses, cal- 
culated mass loss rates peak at max M ~ 0. 1 M^/Gyr, which is just 
large enough to produce the observed transit depths of order 1%. 
Smaller planets with weaker gravities yield larger mass loss rates 
and are also compatible with the observations. 
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Figure 12. Trajectories of the wind in pressure-temperature space, com- 
puted using the full model for M = 0.01, 0.03, and O.OVM®. Gas pressures 
generally remain above the equilibrium vapor pressures of pyroxene, per- 
mitting the condensation of pyroxene grains. 



(ii) Maximum formation mass 

For an assumed planet age of ~5 Gyr, the maximum mass of KIC 
1255b at formation ranges from 0.06 to 0. 1 M^, again depend- 
ing on the energetics of the wind (Figure|9}- 

(iii) Mass threshold for catastrophic evaporation 

A pure rock planet of mass > 0.1 M0 and surface temperature < 
2200 K will survive with negligible mass loss for tens of Gyrs. 

(iv) Time variability' 

The observed occultations of KIC 12557548 vary by up to an order 
of magnitude in depth without any apparent correlation between 
orbits. The implied order-of-magnitude variations in M can be ex- 
plained in principle by the exponential sensitivity of M to condi- 
tions at the planet surface. The dynamical "refresh" timescale of 
the wind ?dyii ~ 14 h is similar to the orbital period Parh = 15.7 h. 
This supports our dusty wind model because were t^y„ » PoA or 
fdyn <K /'orb, then ecllpse depths would correlate from orbit to orbit, 
in violation of the observations. 

(v) Progenitors and occurrence rates 

KIC 1255b's current catastrophic mass loss phase may represent 
only the final few percent of the planet's life. As such, for every 
KIC 1255b-like object, there could be anywhere from 10 to 100 
larger planets in earlier stages of mass loss. These close-in, rela- 
tively quiescent progenitors may be detectable by Kepler through 
conventional "hard-sphere" transits if they are as large as Mercury. 
We cannot, however, rule out the possibility that the progenitors are 
lunar-sized or even smaller. If KIC 1255b remains the only catas- 
trophically evaporating planet in the Kepler database, then the oc- 
cuiTence rate of close-in progenitors orbiting K-stars with sub-day 
periods is > 0.1%, with larger occurrence rates for progenitors in- 
creasingly smaller than Mercury, 
(vi) Iron planet 

KIC 1255b may have lost ~70% of its formation mass to its thermal 
wind. It seems possible that today only the iron core of KIC 1255b 
remains and is evaporating. Transmission spectra of the occulting 



dust cloud might reveal whether the planet's surface is composed 
primarily of iron or silicates, 
(vii) Future modelling and observations 

Keeping the wind hot as it lifts off the planet surface significantly 
enhances mass loss rates. In our model we have included heating 
of the gas by micron-sized grains embedded in the flow. But these 
grains can also shadow the surface from starlight and reduce M 
if they are too abundant. Additional heating may be provided by 
super-blackbody grains with sizes <k 1 pm which do not signif- 
icantly attenuate starlight. Future models should incorporate such 
tiny condensates — and treat the dusty comet-like tail that our pa- 
per has completely ignored. More information about grain size dis- 
tributions and compositions might be revealed by Hubble observa- 
tions of KIC 12557548 scheduled for early 2013. 
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